Having protein alternatives is important for sustainable scaling of agriculture1 to meet the growing food demand. Given the amount of coastal area in many countries, understanding ways to use those resources is important. Marine aquaculture has the potential to play a key role in increasing food production, boosting economic growth in coastal areas, and help maintain clean waterways2.
Marine aquaculture potential is dictated by many variables including ship traffic, dissolved oxygen, bottom depth3. Understanding what coastal areas have potential for different species can guide the implementation of marine aquaculture.
The Exclusive Economic Zones (EEZ) represent areas from coastlines
that countries have resource jurisdiction over, and could target for
developing aquaculture. In this project, I am determing the locations in
West Coast EEZ that are suitable for a variety of oyster aquaculture.
The areas need to fit these specific conditions:
The stars and terra packages will be useful
for handling the raster data.
This analysis will characterize sea surface temperature within the regions with sea surface temperature (SST) from 2008 to 2012. The data is originally generated from NOAA’s 5km Daily Global Satellite Sea Surface Temperature Anomaly v3.1.
This data is stored in the sst_tifs folder: -
average_annual_sst_2008.tif
- average_annual_sst_2009.tif
- average_annual_sst_2010.tif
- average_annual_sst_2011.tif
- average_annual_sst_2012.tif
To characterize the depth of the ocean we will use the General Bathymetric Chart of the Oceans (GEBCO).4
depth.tif for depth dataData on the Exclusive Economic Zones off the west coast of the US are from Marineregions.org.
wc_regions_clean.shp file to access the EEZ
shapefiles# load libraries
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(here)
## here() starts at /Users/lunacatalan/Documents/dev/eds223/final_repos/oyster-aquaculture-suitability
library(stars) # for .tif files
## Loading required package: abind
library(terra)
## terra 1.7.55
##
## Attaching package: 'terra'
##
## The following object is masked from 'package:tidyr':
##
## extract
library(leaflet)
library(patchwork)
##
## Attaching package: 'patchwork'
##
## The following object is masked from 'package:terra':
##
## area
library(shiny)
library(htmltools)
Read in the EEZ and depth data:
# load in data for west coast regions
wc_eez <- st_read(here("data/wc_regions_clean.shp")) %>%
rename(ID = rgn_id)
## Reading layer `wc_regions_clean' from data source
## `/Users/lunacatalan/Documents/dev/eds223/final_repos/oyster-aquaculture-suitability/data/wc_regions_clean.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -129.1635 ymin: 30.542 xmax: -117.097 ymax: 49.00031
## Geodetic CRS: WGS 84
# read in bathymetry ocean depth data
depth <- read_stars(here("data/depth.tif")) %>%
rast()
Since there are multiple .tifs for the sea surface temperature data, I made a list and made a raster stack instead of individually reading in the files and then stacking them. Simplified!
# read in annual average sea surface temperatures
file_names <- list.files("data/sst_tifs/",
full.names = TRUE)
# stack the rasters together
sst_stack <- rast(file_names)
summary(sst_stack)
## average_annual_sst_2008 average_annual_sst_2009 average_annual_sst_2010
## Min. :280.4 Min. :279.4 Min. :280.2
## 1st Qu.:285.3 1st Qu.:285.7 1st Qu.:285.7
## Median :287.1 Median :287.7 Median :287.4
## Mean :287.1 Mean :287.7 Mean :287.5
## 3rd Qu.:289.1 3rd Qu.:289.7 3rd Qu.:289.4
## Max. :301.4 Max. :301.5 Max. :301.0
## NA's :42125 NA's :42204 NA's :42229
## average_annual_sst_2011 average_annual_sst_2012
## Min. :278.9 Min. :278.5
## 1st Qu.:285.5 1st Qu.:285.6
## Median :287.0 Median :287.0
## Mean :287.1 Mean :287.2
## 3rd Qu.:288.8 3rd Qu.:289.0
## Max. :307.3 Max. :310.2
## NA's :42071 NA's :42034
class(sst_stack)
## [1] "SpatRaster"
## attr(,"package")
## [1] "terra"
Check the crs of the data before continuing so that we dont run into trouble when we try to make them interact:
# do the crs match?
st_crs(wc_eez) == st_crs(sst_stack) # no they do not
st_crs(wc_eez) == st_crs(depth) # yes
The crs for the West Coast data and the Sea Surface temp rasters don’t match. We have to re-project so that we can work on them together.
To make sure that we can plot the sea surface temperature with the rest of the data we reproject to the crs of the west coast eez raster.
# reproject sst_stack to same crs as wc
sst_stack <- project(sst_stack,
crs(wc_eez))
# check if the new object has new crs
summary(sst_stack)
## average_annual_sst_2008 average_annual_sst_2009 average_annual_sst_2010
## Min. :280.4 Min. :279.4 Min. :280.2
## 1st Qu.:285.3 1st Qu.:285.7 1st Qu.:285.7
## Median :287.1 Median :287.7 Median :287.4
## Mean :287.1 Mean :287.7 Mean :287.5
## 3rd Qu.:289.1 3rd Qu.:289.7 3rd Qu.:289.4
## Max. :301.4 Max. :301.4 Max. :300.9
## NA's :42125 NA's :42206 NA's :42229
## average_annual_sst_2011 average_annual_sst_2012
## Min. :278.9 Min. :278.5
## 1st Qu.:285.5 1st Qu.:285.6
## Median :287.0 Median :287.0
## Mean :287.1 Mean :287.2
## 3rd Qu.:288.8 3rd Qu.:289.0
## Max. :307.3 Max. :309.9
## NA's :42071 NA's :42034
To use them together we need to make sure that the resolutions and extents of the objects are the same.
The temperature data is in Kelvin. To avoid having to transform the inputs of temperature parameters into Kelvin, we can change the means to Celsius bu substracting by 273.15 degrees.
# mean from stack
sst_mean_K <- terra::mean(sst_stack,
na.rm = TRUE)
# convert values into C by subtracting by 273.15
sst_mean <- sst_mean_K - 273.15
sst_mean
## class : SpatRaster
## dimensions : 480, 408, 1 (nrow, ncol, nlyr)
## resolution : 0.04165905, 0.04165905 (x, y)
## extent : -131.9848, -114.9879, 29.99208, 49.98842 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : mean
## min value : 5.041986
## max value : 32.720087
Use resample and crop from the
terra package to match the resolutions and extents of the
sst and depth data:
# resample depth to match resolution of sst
depth_resample <- resample(depth, sst_mean,
method = "near")
# crop the depth raster to the sst raster
depth_crop <- crop(depth_resample, sst_mean)
Check the resolutions and extents of the depth and sst rasters:
# compare what the extents are visually
plot(depth)
plot(sst_mean)
plot(depth_crop)
# check if the crs match
st_crs(depth_crop) == st_crs(sst_mean)
# look at the resolutions
depth_crop
sst_mean
I am using the sea surface temperature that oysters can suitably grow in: 11-30 deg C - recalssify so that all values outside of the 11-30 degree range are NA
# sea surface temperature: 11-30° C
rcl_temp <- matrix(c(11, 30, 1, # from EQUAL TO 11 to EQUAL TO 30, set the value to 1
-Inf, 11, NA, # from infinity to 11, set the value to 0
30, Inf, NA), # from greater than 30 to infinity, set the value to 0
ncol = 3, byrow = TRUE)
# reclassify
rcl_sst <- classify(sst_mean,
rcl = rcl_temp)
I am using the depth range that oysters can suitably grow of 0-70 meters below sea level. This translates to -70 and 0 meters. Any values outside of this range will be NA.
# depth: 0-70 meters below sea level
rcl_depth <- matrix(c(-70, 0, 1, # from EQUAL TO -70 to EQUAL TO 0, set the value to 1
-Inf, -70, NA, # from infinity to less than -70, set the value to 0
0, Inf, NA), # from greater than 0 to infinity, set the value to 0
ncol = 3, byrow = TRUE)
# reclassify
rcl_depth <- classify(depth_crop,
rcl = rcl_depth)
Since the values in each raster are not the values of 1 (if suitable) and NA (if not suitable), we can perform map algebra and multiply the layers. Anywhere that there are NA’s will return an NA value, and only the pixels that have a value of 1 in both layers will return a 1. Therefore, the resulting raster will only have 1’s in the areas that match the suitability parameters in the depth AND the temp conditions.
# create multiplying function : layers that have 0 will multiply to 0 - only the places where it equals 1 in both will be 1 after multiplying
rast_mult <- function(layer1, layer2) {
(layer1*layer2)
}
#create stack of the reclassified data
depth_sst_stack <- c(rcl_sst, rcl_depth)
# apply the function to each layer in the stack
depth_sst <- lapp(depth_sst_stack, # stacked data
fun = rast_mult) # function
#rename layer as suitable
names(depth_sst) <- "suitable"
# check if it got renamed
depth_sst
## class : SpatRaster
## dimensions : 480, 408, 1 (nrow, ncol, nlyr)
## resolution : 0.04165905, 0.04165905 (x, y)
## extent : -131.9848, -114.9879, 29.99208, 49.98842 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : suitable
## min value : 1
## max value : 1
To work with the EEZ data we can rasterize it using the
raster we created above with suitable areas.
This allows us to mask the EEZ data with the depth_sst
data to identify what areas within the EEZ are suitable. To find the
area that these pixels represent we can use the expanse
function that sums the raster cells that are 1’s and not NA’s.
# rasterize the eez areas
id_rast = rasterize(vect(wc_eez), # create a SpatVector to make into a raster
depth_sst,
field = "rgn") # select layer in database
# cells that are in the right spot
eez_suitable <- mask(id_rast, depth_sst)
# Compute the area covered by polygons or for all raster cells that are not NA in km
eez_zonal <- expanse(eez_suitable,
unit = "km", # output is km2 because its an area
byValue = TRUE)
We can then calculate the percent of each EEZ that is suitable for oyster aquaculture, and comapre it to the total area in each EEZ.
# find the percentage of each zone that is suitable
eez <- wc_eez %>%
mutate(area = eez_zonal$area,
percent = (area/area_km2)*100) # find percent
This gives a guide for regions to understand the productivity/economic potential of the area, and quantifying the total area is good for understanding priority for investing in marine oyster aquaculture.
Create basemap:
# basemap
basemap <- leaflet() %>%
addProviderTiles(
"Esri.WorldImagery",
group = "Esri.WorldImagery"
) %>%
# add a layers control
addLayersControl(
baseGroups = c("Esri.WorldImagery")
)
Leaflet map with basemap of percent suitable areas:
pal <- colorNumeric(
c("#E1F5C4", "#B3E0A6FF", "#7DC370FF", "#59A253FF","#368747FF"),
# colors depend on the count variable
domain = eez$percent,
)
map_percent <- basemap %>%
addPolygons(data = eez,
# set the color of the polygon
color = ~pal(percent),
# set the fill opacity
fillOpacity = 0.9
) %>%
# add a legend
addLegend(
data = eez,
pal = pal,
values = ~percent,
position = "bottomleft",
title = "Percent Suitable \nArea:",
opacity = 0.9
) %>%
addControl("Percent Suitable Area for Oysters In West Coast EEZ",
position = "bottomright",
)
map_percent
pal <- colorNumeric(
c("#F4D166FF", "#F8AF50FF", "#F38C30FF", "#EB6C1CFF","#CB4D22FF"),
# colors depend on the count variable
domain = eez$area,
)
map_total <- basemap %>%
addPolygons(data = eez,
# set the color of the polygon
color = ~pal(area),
# set the fill opacity
fillOpacity = 0.9
) %>%
# add a legend
addLegend(
data = eez,
pal = pal,
values = ~area,
position = "bottomleft",
title = "Total Suitable \nArea:",
opacity = 0.9
) %>%
addControl("Total Suitable Area for Oysters In West Coast EEZ",
position = "bottomright",
)
map_total
Re-set the west coast eez zones file:
wc_eez <- st_read(here("data/wc_regions_clean.shp")) %>%
rename(ID = rgn_id)
## Reading layer `wc_regions_clean' from data source
## `/Users/lunacatalan/Documents/dev/eds223/final_repos/oyster-aquaculture-suitability/data/wc_regions_clean.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -129.1635 ymin: 30.542 xmax: -117.097 ymax: 49.00031
## Geodetic CRS: WGS 84
species_maps :Takes in 5 parameters: species, temp_low,
temp_high, depth_low, and
depth_high. - Will produce two visualizations of maps that
show percent rnak and total area rank for EEZ on the west coast of the
US.
species_maps <- function(species, temp_low, temp_high, depth_low, depth_high) {
# sea surface temperature:deg C
rcl_temp <- matrix(c(temp_low, temp_high, 1,
-Inf, temp_low, NA,
temp_high, Inf, NA),
ncol = 3, byrow = TRUE)
# reclassify
rcl_sst <- classify(sst_mean,
rcl = rcl_temp)
# depth meters below sea level
rcl_depth <- matrix(c(depth_low, depth_high, 1,
-Inf, depth_low, NA,
depth_high, Inf, NA),
ncol = 3, byrow = TRUE)
# reclassify
rcl_depth <- classify(depth_crop,
rcl = rcl_depth)
rast_mult <- function(layer1, layer2) {return(layer1*layer2)}
#create stack of the reclassified data
depth_sst_stack <- c(rcl_sst, rcl_depth)
# apply the function to each layer in the stack
depth_sst <- lapp(depth_sst_stack, # stacked data
fun = rast_mult)
# rasterize the eez areas
id_rast = rasterize(vect(wc_eez), # create a SpatVector to make into a raster
depth_sst,
field = "rgn") # select layer in database
# cells that are in the right spot
eez_suitable <- mask(id_rast, depth_sst)
# Compute the area covered by polygons or for all raster cells that are not NA in km
eez_zonal <- expanse(eez_suitable,
unit = "km", # output is km2 because its an area
byValue = TRUE) %>%
rename(rgn = value)
# find the percentage of each zone that is suitable
eez <- wc_eez %>%
right_join(eez_zonal, by = "rgn", copy = TRUE) %>%
mutate(percent_area = (area/area_km2)*100) # find percent
pal <- colorNumeric(
c("#E1F5C4", "#B3E0A6FF", "#7DC370FF", "#59A253FF","#368747FF"),
# colors depend on the count variable
domain = eez$percent,
)
percent <- leaflet(eez) %>%
addPolygons(data = eez,
# set the color of the polygon
color = ~pal(percent_area),
# set the fill opacity
fillOpacity = 0.9) %>%
addLegend(
data = eez,
pal = pal,
values = ~area,
position = "bottomleft",
title = "Total Suitable \nArea:",
opacity = 0.9) %>%
addControl(paste("% Suitable Area <br> for", species), # updates with species input,
position = "bottomright",
) %>%
addScaleBar() %>% # add scalebar
addTiles() # add basemap
pal <- colorNumeric(
c("#F4D166FF", "#F8AF50FF", "#F38C30FF", "#EB6C1CFF","#CB4D22FF"),
# colors depend on the count variable
domain = eez$area,
)
map_total <- basemap %>%
addPolygons(data = eez,
# set the color of the polygon
color = ~pal(area),
# set the fill opacity
fillOpacity = 0.9
) %>%
# add a legend
addLegend(
data = eez,
pal = pal,
values = ~area,
position = "bottomleft",
title = "Total Suitable \nArea:",
opacity = 0.9
) %>%
addControl(paste("Total Suitable Area <br> for", species), # title that updates with species input
position = "bottomright") %>%
addScaleBar() %>%
addTiles()
# put the leaflet maps next to each other
map_grid <-
tagList(
tags$table(width = "100%",
tags$tr(
tags$td(percent),
tags$td(map_total)
)
)
)
# show the maps
return(map_grid)
}
# testing the function
species_maps("Clam", 6, 12, -200, -30)
Fisheries, N. (2022, September 15). Aquaculture supports a sustainable earth. https://www.fisheries.noaa.gov/feature-story/aquaculture-supports-sustainable-earth↩︎
Aragão, Cláudia, et al. “Alternative Proteins for Fish Diets: Implications beyond Growth.” Animals : An Open Access Journal from MDPI, U.S. National Library of Medicine, 7 May 2022, www.ncbi.nlm.nih.gov/pmc/articles/PMC9103129/.↩︎
GEBCO Compilation Group (2022) GEBCO_2022 Grid (doi:10.5285/e0f0bb80-ab44-2739-e053-6c86abc0289c).↩︎
GEBCO Compilation Group (2022) GEBCO_2022 Grid (doi:10.5285/e0f0bb80-ab44-2739-e053-6c86abc0289c).↩︎